import random
import numpy
#import matplotlib
#from matplotlib import pyplot as plt
import isotopes
import math
import time
import pickle

verbose = True #Prints out all kinds of debug crap

interval_time = 3*3600   #seconds
intervals = 100          # num of intervals to save
chns = 1024              # chanels in MCA
energy_scale = 80.       #chn/MeV
energy_offset = 0.1      #MeV @ chn=1

# dead layer thickness
h = 0.5 #um
dEdX = 0.11 #MeV/um
gauss_sigma = 0.009 #MeV


ESC_volume = 10./1000.  #m^3
LOOP_volume = 2./1000.  #m^3

# Initial conditions (typical?)
emmination_rate = {}       #Sample  #Chamber  #Loop
emmination_rate['Rn222'] = (50.,    1.,       1.)
#emmination_rate['Rn220'] = (0.0,    0.,       1.)
#emmination_rate['Rn219'] = (0.0,    0.,       1.)

# diode_rate_Pb210 = 100. #dpd
t0 = {}
t0['Rn222']   = 100    #atoms
#t0['Rn220']   = 0    #atoms
#t0['Rn219']   = 0    #atoms

# This is very slow way to get Pb210 events on the diode, I dont recommend using it
t0_diode = {}
t0_diode['Pb210']  = 0  #atoms
# Non probabilistic way to add Pb210 events
Pb210_diode_dpd = 100

col_eff = {}
col_eff['Po218'] = 0.8
col_eff['Pb214'] = 0.8

def main():
   # Setup counters for isotopes and initial conditions
   ESC_contents = {}
   for iso in t0.keys():
      ESC_contents[iso] = t0[iso]
   DIODE_contents = {}
   for iso in t0_diode.keys():
      DIODE_contents[iso] = t0_diode[iso]
   
   run_events = []
   
   # Start running MC, @ resolution of 1 second
   for interval in range(intervals):

      diode_inject = []
      if verbose: print "Interval=%d"%interval

      for dt in range(interval_time):
      
         # How many emminations occur based on emmination
         for iso in emmination_rate.keys():
            if iso not in ESC_contents.keys():
               ESC_contents[iso] = 0
            ESC_contents[iso] += sum(numpy.random.random(sum(emmination_rate[iso]))<(1./3600./24.))
         
         # Decay atoms on the DIODE
         for iso in DIODE_contents.keys():
            iso_decays = int(sum(numpy.random.random(DIODE_contents[iso]).astype('float128') <
               numpy.longdouble(1.-numpy.exp(-0.693147/isotopes.decays[iso][0]['T12']))))

            # Subtrace decayed parents
            DIODE_contents[iso] -= iso_decays
            
            # determine which branch each decaying atom takes
            for id in range(iso_decays):
               rand_br = numpy.random.random()
               tbr = 0.
               for branch in isotopes.decays[iso]:
                  tbr+= branch['BR']
                  if tbr > 1.0: print "Branching ratio is >1.0!!"
                  if rand_br < tbr:
                     if branch['daughter'] == None:
                        break 
                     # It takes this branch
                     if branch['daughter'] not in DIODE_contents.keys():
                        # Check if that atom exists yet in the ESC, If not create the dict entry
                        DIODE_contents[branch['daughter']] = 0
                     DIODE_contents[branch['daughter']] += 1
                
               # And determine the energy that a alpha imparts into the diode
               if branch['type'] == 'alpha':
                  diode_inject.append(diode_responce(branch['energy']))
         
         
         # How many atoms decay in the ESC + loop?
         # Loop over all the atoms we have presently and see if they decay
         for iso in ESC_contents.keys():
            
            iso_decays = int(sum(numpy.random.random(ESC_contents[iso]).astype('float128') <
               numpy.longdouble(1.-numpy.exp(-0.693147/isotopes.decays[iso][0]['T12']))))
            
            # Subtract decayed parents
            ESC_contents[iso] -= iso_decays
            
            # determine which branch each decaying atom takes
            for id in range(iso_decays):
               rand_br = numpy.random.random()
               tbr = 0.
               for branch in isotopes.decays[iso]:
                  tbr+= branch['BR']
                  if tbr > 1.0: print "Branching ratio is >1.0!!"
                  if rand_br < tbr:
                     # It takes this branch
                     if branch['daughter'] == None:
                        break
                     if branch['daughter'] not in ESC_contents.keys():
                        # Check if that atom exists yet in the ESC, If not create the dict entry
                        ESC_contents[branch['daughter']] = 0
                     ESC_contents[branch['daughter']] += 1
                     break

               # Determine if the daughter will drift to the diode
               if branch['daughter'] in col_eff.keys() and branch['daughter'] != None:
                  rand_col = numpy.random.random()
                  if rand_col < col_eff[branch['daughter']]:
                     if iso not in DIODE_contents.keys():
                        DIODE_contents[iso] = 0
                     DIODE_contents[iso]+=1
                     ESC_contents[branch['daughter']] -= 1
      
      # Inject a constant amount of Bi210 into the diode as well
      for i in range(int(math.ceil(Pb210_diode_dpd*interval_time/(24.*3600.)))):
         diode_inject.append(diode_responce(isotopes.decays['Po210'][0]['energy']))
      print interval, ESC_contents
      print diode_inject
   
      run_events.append(diode_inject)

   # Output data in dictionary
   ESCSim = {}
   ESCSim['events'] = run_events[:]
   ESCSim['emmination_rate'] = emmination_rate
   ESCSim['interval_time'] = interval_time
   ESCSim['intervals'] = intervals
   ESCSim['chns'] = chns              # chanels in MCA
   ESCSim['energy_scale'] = energy_scale       #chn/MeV
   ESCSim['energy_offset'] = energy_offset      #MeV @ chn=1

   ESCSim['ESC_volume']    = ESC_volume  #m^3
   ESCSim['LOOP_volume']   = LOOP_volume  #m^3

   ESCSim['t0']         = t0
   ESCSim['t0_diode']   = t0_diode

   ESCSim['col_eff'] = col_eff
   
   pickle.dump(run_events,open('ESCSim_%s.pkl'%time.strftime('%m%d_%H%M',time.localtime()),'w'))
               
def diode_responce(e_in):
   if e_in == None:
      return 0
   
   theta = numpy.random.random() # 2*pi
   if theta >= 0.5:
      #skip it, alpha went the wrong way
      return 0
   phi = math.acos(2.*numpy.random.random() - 1.)

   #The length the track is in the silicone
   L = e_in/dEdX

   # the dead layer wont make detectable energy
   L_dead = h/math.sin(phi)

   #Energy remaining to deposit in active layer
   E_dep = (L-L_dead)*dEdX
   #if E_dep < 5750:
   #   continue
   
   #Smear with a gaussian
   E_dep_smear = numpy.random.normal(E_dep,gauss_sigma,1)[0]
   #chns = 1024              # chanels in MCA
   #energy_scale = 80.       #chn/MeV
   #energy_offset = 0.1      #MeV @ chn=1
   
   return int(math.floor((energy_offset+E_dep_smear)*80.))
               

if __name__ == '__main__':
   main()
